*___________________________________________________________________________________________________

	VarianceDecomp 1.00
	
	Joel Hasbrouck
	Department of Finance
	Stern School of Business
	New York University
	Mail Code 0268
	44 West 4th St.
	New York, NY 10012

	email: jhasbrou@stern.nyu.edu

	Computes a variance decomposition. You may tell it to use a particular permutation,
	but it does not examine all decompositions. 

	Keyword Parameters:

	dsCov		Input dataset that contains the covariance matrix of disturbances.
				The only numeric variables in this dataset should be the covariances.
				(If there are others, drop them before calling this routine.)
	dsCoeff		Input dataset that contains the coefficient matrix.
				A decomposition will be computed for each row of the matrix.

	Optional parameters
	RowLabel	Name of input dataset variable that names the rows of the coefficient matrix.
	Perm		Permutation of variables to use.
	PrintLevel	1	Summary printed output (default)
				2	Detailed output

	A test driver program follows the macro definition.
____________________________________________________________________________________________________;
%macro VarianceDecomp(dsCov=, dsCoeff=, RowLabel=, Perm=, PrintLevel=2);
%*	Get number of variables;
title 'Variance Decomposition';
data _null_;
	if 0 then set &dsCov nobs=countxx;
	call symput("n",compress(put(countxx,12.)));
	stop;
	run;
proc iml;
	start main;
	reset printadv=1;
	use &dsCoeff;
	%if %length(&RowLabel)=0 %then %do;
		read all into b [colname=VarNamesb];
		RowLabel = concat('Coefficient vector', char(1:nrow(b),2,0));
	%end;
	%else %do;
		read all into b [colname=VarNamesb rowname=&RowLabel];
		RowLabel = &RowLabel;
	%end;
	print b [rowname=RowLabel label="Coefficient matrix"];
	n = ncol(b);
	nRowb = nrow(b);
	use &dsCov;
	read all into cov [colname=VarNames];
	print cov [rowname=VarNames label="Covariance matrix"];
	if n^=nrow(cov) then do;
		print "Covariance matrix not square or not conformable with coefficients.";
		abort;
	end;

	sd = diag(vecdiag(cov)##(-1/2));
	cor = sd * cov * sd;
	print cor [rowname=VarNames colname=VarNames label='Correlation matrix' format=6.3];

	*	Set vector of variable names that describes the permutation;
	perm=VarNames%str(;);	*(default);
	%if %length(&perm)>0 %then %do;
		%do i=1 %to &n; 
			perm[&i] = "%scan(&perm,&i)"%str(;)
		%end;
		if ncol(perm)^=n then do;
			print "Supplied permutation vector does not include all variables:";
			print perm;
			abort;
		end;
	%end;

	*	Set vector of indexes for permutation;
	iperm = 1:n;
	do i=1 to n;
		iperm[i]=loc(perm[i]=varNames);
	end;
	varNamesP = PermuteCols(varNames,iperm,1);
	print VarNamesP [label="Purmutation used in decomposition / ordering of variables:"];

	*	Execute the permutations;
	bP = PermuteCols(b,iperm,1);
	%if &PrintLevel>=2 %then
		print bp [rowname=RowLabel colname=varNamesP label="Permuted coefficients"];;
	covP = PermuteMatrix(cov,iperm,1);
	%if &PrintLevel>=2 %then
		print covP [rowname=varNamesP colname=varNamesP label="Permuted covariance matrix"];;

	*	Analyze the variance contributions;
	cd = t(root(covP));
	print cd [label='Cholesky factor of permuted covariance matrix' rowname=RowLabel colname=VarNamesP];
	VarContrib = (bP*t(root(covP)))##2;
	print VarContrib [rowname=RowLabel colname=varNamesP label="Variance contributions (ordered)"];
	VarTotal = VarContrib[,+];
	print VarTotal [rowname=RowLabel label="Total variance"];
	VarTotal = t( shape(VarTotal,ncol(VarContrib), nrow(VarContrib) ));
	VarProp = VarContrib/VarTotal;
	print VarProp [rowname=RowLabel colname=varNamesP label="Proportional contributions" format=6.3];
	finish main;

	*______________________________________________________________________________________________;
	start PermuteCols(a,p,direction);
	*	Reorders the columns of a matrix.
		direction=+1 for the forward permutation and -1 to reverse the permutation;
	aP = a;
	do i=1 to ncol(a);
		k = p[i];
		if direction=1 then aP[,i] = a[,k];
		else aP[,k] = a[,i];
	end;
	return (aP);
	finish PermuteCols;

	*______________________________________________________________________________________________;
	start PermuteMatrix(a,p,direction);
	*	Reorders the variables in preparation for Cholesky decomposition.
		direction=+1 for the forward permutation and -1 to reverse the permutation;
	aP = a;
	if nrow(a)=1 | ncol(a)=1 then do i=1 to max(ncol(a),nrow(a));
		k = p[i];
		if direction=1 then aP[i] = a[k];
		else aP[k] = a[i];
	end;
	else do i=1 to nrow(a);
		do j=1 to ncol(a);
			k = p[i];
			l = p[j];
			if direction=1 then aP[i,j] = a[k,l];
			else aP[k,l] = a[i,j];
		end;
	end;
	return (aP);
	finish PermuteMatrix;

	run;
	quit;
%mend VarianceDecomp;

*___________________________________________________________________________________________________

	Test program:
____________________________________________________________________________________________________;
/*options nocenter mprint;*/
/**	Generate some correlated variates;*/
/*title 'Decomposition with simulated data';*/
/*data;*/
/*	do i=1 to 100;*/
/*		p1 = normal(123);*/
/*		p2 = p1 + normal(456);*/
/*		p3 = p2 - normal(789);*/
/*		p4 = p2 - normal(246);*/
/*		p5 = p3 + normal(135);*/
/*		output;*/
/*		drop i;*/
/*	end;*/
/*proc corr cov outp=acov;*/
/*data acov;*/
/*	set acov;*/
/*	where _type_="COV";*/
/*proc print data=acov;*/
/*data coeff;*/
/*	p1 = .4;*/
/*	p2 = .8;*/
/*	p3 = .2;*/
/*	p4 = -.1;*/
/*	p5 = .5;*/
/*	type = "Coefficients A";*/
/*	output;*/
/*	p5 = .25;*/
/*	type = "Coefficients B";*/
/*	output;*/
/*	run;*/
/*%VarianceDecomp(dsCov=acov, dsCoeff=coeff, RowLabel=type, perm=p2 p3 p1 p5 p4, PrintLevel=1);*/
/*%VarianceDecomp(dsCov=acov, dsCoeff=coeff, RowLabel=type, perm=p2 p3 p1 p5 p4, PrintLevel=1);*/
